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INTRODUCTION  • 


Andrew  and  Fleming  (ref  1)  discussed  the  calculation  of  null  geodesics  for  the 
Schwarzchild,  Kerr-Newman,  and  Winicour  metrics.  Numerical  computation  employed 
FORTRAN  code,  which  was  produced  by  Mathematica  (ref  2).  The  Mathematica  code  employed 
to  generate  the  FORTRAN  code  for  the  neutral  (Q  =  0)  Kerr-Newman  geodesics  was  presented. 
Reference  1  provides  a  convincing  demonstration  of  the  utility  of  symbolic  manipulation  software 
for  the  development  of  error-free,  complex  FORTRAN  (or  C)  code. 

In  this  report,  a  general  approach  to  problems  in  tensor  analysis  (ref  3)  employing 
Mathematica  is  described.  The  standard  expressions  of  tensor  calculus  are  transcribed  directly 
into  Mathematica  modules  in  Part  I.  The  operation  of  the  modules  is  illustrated  as  they  are 
described  by  application  to  a  simple  two-dimensional  curved  space.  The  present  formulation 
allows  enumeration  and  simplification  of  complex  tensor  equations,  employing  built-in 
Mathematica  functions  (such  as  Expand[],  Together[],  or  Simplify[]). 

The  code  is  applied  to  address  some  simple  problems  pertinent  to  Schwarzchild  space- 
time  in  Part  II.  Employing  the  built-in  Mathematica  operator  Simplify[],  the  modules  described 
produced  FORTRAN  coded  Runge-Kutta  geodesic  equations  for  Kerr-Newman  space-time  in 
less  than  thirty  minutes  on  a  386-based  PC  operating  at  30  MHz.  Of  course,  if  simplification 
based  on  more  efficient  Mathematica  operations  (e.g.,  combinations  of  Apart[],  Expand[],  and 
Together[])  can  be  achieved,  symbolic  computation  times  can  be  dramatically  reduced. 
Computation  times  are  much  shorter  for  the  Schwarzchild  geodesic  equations. 


PART  I.  TENSOR  ANALYSIS  IN  Mathematica 
Covariant  Metric  Tensor 


The  starting  point  for  a  typical  problem  in  tensor  analysis  is  the  specification  of  the 
(covariant)  metric  tensor. 


The  operation  of  the  Mathematica  tensor  analysis  modules  defined  here  is  demonstrated 
by  applications  to  geometry  on  the  surface  of  a  helicoid  in  R^  having 


g  =  = 


(\ 

lo 


0  ) 


where  the  coordinates  are  x[l][s]  and  x[2][s]  and  c  is  constant.  This  g  is  referred  to  as  the 
"example  metric"  in  Part  I.  One  encodes  this  g  directly  into  Mathematica  via 


g=DiagonalMatrix[{l,c^2  -I-  x[l][s]^2}]; 

Note  that  coordinates  are  given  in  the  form  x[i][s].  For  the  current  analysis,  one  could  simply 
use  x[i].  However,  as  in  Part  II,  one  is  often  interested  in  derivatives  with  respect  to  the  line 
element  s,  and  we  choose  to  include  this  dependence  from  the  beginning. 
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Christoffel  Symbols 


The  first  step  in  a  tensor  problem  is  the  computation  of  the  Christoffel  symbols.  ’ 
Christoffel  symbols  of  the  first  kind  are  defined  as 


^8a 


+ 


ac“ 


dx‘  j 


Christoffel  symbols  of  the  second  kind  are  given  by 

lip  =  Es*‘r,.p  =  s‘T,.p 

i 


where  the  last  equality  introduces  "summation  convention"  which  is  adopted  in  this  report  and 
the  contravariant  metric  tensor  is  the  inverse  of  the  covariant  metric,  i.e., 

The  following  Mathematica  module  defines  the  contravariant  metric  tensor  and 
Christoffel  symbols  for  arbitrary  metric  tensors  in  spaces  of  arbitrary  dimension. 

setup[g_]  ;=  Modul  e  [{  dim = Length  [g] } , 
ginv  =Together[Inverse[g]]; 

GAMMA  =Array[gl,{dim,dim,dim}]; 
gamma  =Array[g2,{dim,dim,dim}]; 

Do[  Do[  Do[GAMMA[[i,k,j]]=GAMMA[[i,j,k]]= 

Together[(l/2)(D[g[[i,j]],x[k][s]]  + 

D[g[[k,i]],x[j][s]]-D[g[[j,k]],x[i][s]])], 

{i,dim}],{k,j,dim}],{j,dim}]; 

Do[  Do[  Do[gamma[[i,k,j]]=gamma[[i,j,k]]= 

Together[Sum[ginv[[i,l]]GAMMA[[l,j,k]],{l,dim}]], 

{i,dim}],{k,j,dim}],{j,dim}]  ] 

N.  B.,  the  partial  derivative  of  a  function  f[u,v,...,x,...]with  respect  to  the  variable  x  is  computed 
via  D[f[u,v,...,x,...],x]in  Mathematica.  N.B.,  the  Mathematica  Together[]  operator  was  effective  in 
simplifying  the  resulting  expressions  for  a  number  of  "simple  problems."  Of  course,  the  user 
could  apply  an  other  built-in  (e.g.,  Simplify[])  or  user-defined  Mathematica  operator. 

After  setup[g]  is  run,  the  contravariant  metric  tensor  is  returned  as  ginv,  Christoffel 
symbols  of  the  first  kind  are  in  the  array  GAMMA,  i.e.,  =  GAMMA[[i,j,k]],  and  Christoffel 

symbols  of  the  second  kind  are  in  the  array  gamma,  i.e.,  F  ‘  =  gamma[[i,j,k]]. 
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To  compute  the  Christoffel  symbols  and  g'^  for  the  test  metric,  enter 
setup[g]. 

Example:  Display  gamma  for  the  example  metric, 
gamma 


returns 

x[l][s]  x[l][s] 

{{{0,  0},  {0,  -x[l][s]}},{{0, - },  { - ,  0}}}, 

2  2  2  2 
c  +x[l][s]  c  +x[l][s] 

{{{0,  0},  {0,-x[l][s]}}, 

Exercise:  Code  a  Mathematica  module  to  compute  covariant  derivatives. 

The  Curvature  Tensor 


The  curvature  tensor  can  be  computed  directly  from  the  Tj^  ^ 

D  •  _  •pP  j4  _  pP  p* 

■  TT  “  T7  ^  “y^  P*  ^  «* 


ax' 


dx^ 


Transcription  of  the  expressions  for  the  is  straightforward, 

curvature[gamma_,i_,a_,j_,k_]: =Module[{b,dim =Length[gamma]}, 
Together[ 

D[gamma[[i,a,j]],x[k][s]]-D[gamma[[i,a,k]],x|j][s]]  + 
Sum[gamma[[b,a,jj]gamma[[i,b,k]]- 

gamma[[b,a,k]]gamma[[i,b,j]],{b,dim}]  ]] 

and  full  curvature  tensor  is  returned  by 


curvature  [gamma_] :  =  Module  [  { dim = Length  [gamma] } , 
B=Array[rl,{dim,dim,dim,dim}]; 
Do[Do[Do[Do[B[[i,a,j,k]]=curvature[gamma,i,a,j,k], 
{i,dim}],{a,dim}],{j,dim}],{k,dim}];B] 


The  covariant  form,  is  called  the  Riemann-Christoffel  curvature  tensor. 
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Example:  Display  the  curvature  tensor  for  the  example  metric 
B  =  curvature[gamma] 


returns 

2  2 
c  c 

{{{{0,  0},  {0,0}},  {{0, - },{-( - ),  0}}}, 

2  2  2  2 

c  +x[l][s]  c  +x[l][s] 


2  ^  2 
c  c 

{{{0,  -( - )},  { - ,  0}},{{0,  0},  {0,  0}}}} 

2  22  2  22 
(c  +x[l][s])  (c  +x[l][s]) 


The  Ricci  Tensor 

The  Ricci  tensor  plays  an  important  role  in  relativity  theory.  It  is  defined  as  a 
contraction  of  the  curvature  tensor. 


,  T'P  p«  - 


Thus,  to  Print}]  the  components  of  R-  =  B°  enter  the  Mathematica  statement 
Do[Do[Print["R(",i,",",j,")  =  ",Simplify[  Sum[B[[a,i,a,j]],{a,2}]  ]  ], 

{i,2}],  {j,2}] 

which  yields 

2 

R(l,l)  =  ....  J.. - 

2  2  2 

(c  +  x[l][s] ) 

R(2,l)  =  0 
R(l,2)  =  0 

2 

c 

R(2,2)  = - 

2  2 

c  +x[l][s] 

N.B.,  Sum[cj^/'C55/ort,{i,n}]  returns  the  sum  of  expression  evaluated  for  i-values  running  from  1 
to  (the  integer)  n.  'Do[expression,{i,n}]  works  similarly. 
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We  also  transcribe  separate  Mathematica  modules  to  compute  the  Ricci  tensor, 

Ricci[gamma_,i_,j_]:=Module[{m,n,dim=Length[gamma]}, 

T0gether[Sum[D[gamma[[m,i,m]],x[j][s]]-D[gamma[[m,i,j]],x[m][s]]+ 

Sum[gamma[[n,i,m]]gamma[[m,n,j]]- 

gamma[[n,i,j]]gamma[[m,n,m]],{n,dim}],{m,dim}]]] 

and 

Ricci  [gammaj:  =Module[{dim = Length[gamma] }, 

R=Array[rl,{dim,dim}]; 

Do[Do[R[[i,j]]=Ricci[gamma,i,j],{i,dim}],{j,dim}];R] 

Example:  Direct  computation  of  Ricci  tensor  for  the  example  metric. 

Ricci  [gamma] 


returns 


2 


2 


c  c 

{{ - ,  0},  {0, - }} 

2  2  2  2  2 
(c  +  x[l][s]  )  c  +  x[l][s] 

which  is  identical  to  the  results  obtained  by  contraction  of  B,  as  expected. 
Exercise:  Compute  the  Christoffel  symbols  and  Ricci  tensor  for 


1  0  0 
0  1  0 
0  0  1 
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PART  II.  APPLICATION:  THE  SCHWARZCHILD  METRIC  . 


Development  of  a  spherically-symmetric  Ijg^- II  consistent  with  the  Einstein  equations. 
Assumption:  ljg',^JI  is  spherically-symmetric.  Assume  that 


g  '  l*«l  = 


m 

0 

0 

0 


0  0  0^ 
0  0 
0  -r^sin^CO)  0 
0  0  M[rl 


0  0  O' 

0  -x[l][sf  0  0 

0  0  -jc[1][5]W(42]M)  0 

,0  0  0  AfI41]M]> 


where  the  functions  L[r]  and  M[r]  are  to  be  determined.  Encode  via,. 
g=DiagonalMatrix[{L[x[l][s]],-x[l][s]  ^  2,-(x[l][s]Sin[x[2][s]])  ^  2,M[x[l][s]]}]; 

Conventional  Notation 


Direct  substitution  is  achieved  in  Mathematica,  via  the  construction 
expression!  .\x->'w 

which  returns  expression  with  occurrences  of  u  replaced  by  v.  Thus,  to  present  output  in 
{r,q,j,t}-form,  define  the  following  replacement  rules: 

rt={x[l][s]->r,x[2][s]->q,x[3][s]->j,x[4][s]->t}; 

N.B.,  Mathematica  2.1  does  not  have  Greek  fonts.  We  substitute  j  for  phi  and  q  for  theta  here 
and  in  the  sequel  with  a  text  editor. 

Compute  the  Christoffel  symbols,  etc. 

setup[g]; 
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The  Einstein  Equations 


The  Einstein  equations  require  that  Ry  =  0.  Compute  and  display  the  Ricci  tensor. 
R=Simplify[Ricci[gamma]]/.rt 

returns  the  Ricci  tensor  for  the  spherically-symmetric  metric  function, 

2 

L’[r]  U[r]  M’[r]  M’[r]  M”[r] 

{{-(- - ) . + - ,  0,  0,  0}, 

rL[r]  4L[r]M[r]  2  2  M[r] 

4M[r] 

1  r  L’[r]  r  M’[r] 

{0,  -1 . + . ,  0,  0}, 

L[r]  2  2  L[r]  M[r] 

2L[r] 

2  2 
{0,  0,  (Sin[T]  (-2  L[r]  M[r]  -  2  L[r]  M[r]  + 

2 

r  M[r]  U[r]  -  r  L[r]  M’[r]))  /  (2  L[r]  M[r]) ,  0}, 

2 

M’[r]  U[r]M’[r]  M’[r]  M”[r] 

{0,  0,  0,  — . — . . + - }} 

r  L[r]  2  4  L[r]  M[r]  2  L[r] 

4L[r] 

Solution  of  the  Einstein  Equations 

We  seek  functions  L[r]  and  M[r]  such  that  the  Einstein  equations  are  satisfied. 

Step  1.  Eliminate  the  M”[r]  terms  from  the  1,1  and  4,4  parts  of  Rjj  =  R[[i,j]]  =  0. 
temp=((Simplify[(R[[l,l]]M[r[s]]-R[[4,4]]L[r[s]])/.rt])==0) 

returns 

M[r]  L’[r]  +  L[r]  M’[r] 

.( - )  ==  0 

rL[r] 
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Step  2.  Solve  for  L[r]  in  terms  of  M[r].  Use 


WgijW  - 


and  the  Mathematica  operator  DSolve[]  to  obtain  RuleO: 

RuleO = ExpandAll  [Flatten[DSolve[ 

{L[Infinity]  =  =-l,M[Infinity]==l,temp},L[r],r]]] 

Mathematica  issues  a  "warning:" 

Solve::ifun: 

Warning:  Inverse  functions  are  being  used  by  Solve,  so  some  solutions  may  not  be  found. 

then  returns  the  solution  of  interest 

1 

|L[r]  ->  -(- - )} 

M[r] 

Step  3.  Mathematica  will  not  apply  the  rule  for  L[r]  to  replace  L’[r].  Thus,  one  specifies 
a  separate  rule  for  L’[r]  and  defines  Rulel  as  follows: 

Rulel=Flatten[{RuleO,Solve[temp,L’[r]]/.RuleO}]/.r->x[l][s] 

which  returns 

1  M’[x[l][s]] 

{L[x[l][s]]  ->  -( - ),  L’[x[l][s]]  -> - } 

M[x[l][s]]  2 

M[x[l][s]] 

N.B.,  Flatten[{{u},{v}}]  returns  {u,v},  etc.  Thus,  the  subrules  comprising  Rulel  will  be 
applied  consecutively. 

Step  4.  Use  Rulel  and  R22  =  0  to  define  the  form  of  M.  Express  M  in  a  rule  (Rule2). 

Rule2=Flatten[DSolve[ 

0==Factor[Expand[R[[2,2]]/.RuIel]]/.rt,M,r]] 

which  returns  the  rule  governing  the  function  M, 


-10  0  0 
0-100 
0  0-10 
0  0  0  1 
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qi] 

{M  ->  (1  +  — -  &  )} 

#1 

where  C[l]  is  a  constant  of  integration.  The  rule  on  the  function  works  as  expected,  i.e., 
M[x[l][s]]/.Rule2 

returns 

qi] 

1  + - 

x[l][s] 

The  Schwarzchild  Metric 

Applying  Rulel  and  Rule2,  one  obtains  the  form  of  the  spherically-symmetric  metric 
tensor  consistent  with  R  =  0: 

gw=g/.Rulel/.Rule2/.rt 

returns 

1  2 
{{.( - ),  0,  0,  0},  {0,  -r  ,  0,  0}, 

qi] 

1  +  — - 


2  2  C[l] 

{0,  0,  -(r  Sin[q] ),  0},  {0,  0,  0, 1  +  --}} 

r 

As  expected,  this  is  the  Schwarzchild  metric  form.  The  usual  form  of  the  metric  tensor  is 
obtained  by  identifying  the  integration  constant  C[l]  with  -2  G  M/c^  where  G  is  the  gravitational 
constant,  M  is  the  central  mass,  and  c  is  the  speed  of  light. 

Exercise:  Demonstrate  that  R  =  0  for  Schwarzchild  gravitation. 

Exercise:  The  geodesic  trajectories  are  given  by  solutions  of  the  geodesic  equations 

dh[i][s]  ,  w  ^[a]M  dxmis]  _  n 
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Derive  the  geodesic  equations  for 


1  0  0 
llgyll  =010 
0  0  1 

Exercise:  Derive  the  geodesic  equations  for  the  "example  metric," 

^  ®|| 

0  :c[l][jf +c2|| 

Exercise:  Planetary  orbits. 

1.  Derive  the  geodesic  equations  for  the  Schwarzchild  metric. 

2.  Demonstrate  that  if  q[0]  =  p/2  and  q’[0]  =  0,  then  q[s]  =  p/2. 

a.  Express  the  geodesic  equations  for  q[s]  =  p/2. 

3.  Show  that  r[s]  ^  2  j"[s]  =  h  =  constant. 

4.  Show  that  (1  -  (2G  M/c^)  /  r[s])  t’[s]  =  K  =  constant. 

5.  Use  the  results  of  parts  1  through  4  and 

1  _  dx\jm 

ds^  ds  ds 

4 

=  8ijxm'[s]x\jr[s}  -  Si/cm'[sf 

i=l 

to  derive  the  Mathematica  expression 


2  2  2 
K  h  (x[l])’[s] 


C[l]  2  C[l] 

1  + -  x[l][s]  1  + - 

x[l][s]  x[l][s] 


6.  Make  the  change  of  variables, 
x[l][s]->l/u[x[3][s]], 

and  employ  the  condition  derived  in  3,  to  establish  that 

(x[l])’[s]  ->  -(h  u’[x[3][s]]) 
and  derive  the  equations: 

2  2  2 
2  2  K  h  u’[phi] 

1  ==  -(h  u[phi] )  + - 

1  +  C[l]  u[phi]  1  +  C[l]  u[phi] 

2 

2  2  -2  K  C[l]u[phi]  3 

u’[phi]  +  u[phi]  ==  -h  + . C[l]  u[phi] 

2  2 
h  h 


C[l]  3C[l]u[phi] 

u”[phi]  +  u”[phi]  ==  . 

2  2 
2h 

The  last  equation  for  u[phi]  (=  l/r[j])  frequently  serves  as  the  basis  for  a  discussion  of  the 
advance  of  the  perihelion  of  planetary  orbits.  See,  for  example,  Crandall  (ref  4).  How  should 
this  approach  be  modified  to  treat  photon  trajectories? 

Exercise:  Compute  the  Christoffel  symbols,  curvature  tensor,  Ricci  tensor,  and  geodesic 
equations  for  charge-free  (i.e.,  Q  =  0)  Kerr-Newman  gravitation.  Demonstrate  that  R  =  0  for 
the  Kerr-Newman  metric. 
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